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1.   INTRODUCTION 

There  is  no  doubt  that  the  available  technological  possibilities  could, 
justify  hardware  implementation  of  a  much  wider  class  of  functions  than  is 
presently  done.   If  the  corresponding  algorithms  have  similarity,  this  is 
even  more  true.   One  effective  way  to  obtain  such  a  class  of  algorithms  is  to 
use,  in  a  convenient  way,  continued  products  (CP)  or  continued  sums  (CS)  during 
the  function  evaluation.   The  use  of  continued  products  in  the  calculation  of 
some  elementary  functions  appears  as  early  as  1959*  in  Voider' s  CORDIC 
technique[  5] •   The  main  results  of  this  approach  have  been  recently  summarized 
in  the  form  of  a  unified  algorithm  by  ¥alther[6].   Specker[7]  derived  also  a 
class  of  algorithms  using  the  concept  of  continued  products.   Of  particular 
importance  and  usefulness  are  the  results  obtained  by  DeLugish[  1] . ■  He  has 
defined  efficient  algorithms  for  a  wide  class  of  functions  including  division, 
multiplication,  square  root,  logarithm,  exponential,  trigonometric  and  inverse 
trigonometric  functions,  with  operation  times  from  1  to  3  multiplication 
cycle  times.   These  algorithms  are  specified  for  the  radix  2,  using  a  redundant 
digit  set  [-1,  0,  1}  in  continued  products  (sums).   The  main  idea  is  replacement 
of  a  required  operation  or  function  evaluation  by  two  simple  step  by  step 
processes  using  addition/subtraction,  shifting  and  possibly  a  set  of  precomputed 
constants,  stored  in  a  read-only  memory.   One  of  the  processes  is  normalization, 
through  which  the  digits  of  continued  products  (sums)  are  generated  and  another 
is  the  related  result  evaluation.   These  processes  can  be  carried  out  in 
parallel,  so  for  a  fast  operation  at  least  two  arithmetic  units,  almost  the 
same,  are  required. 


The  work  done  here  is  based  upon  the  results  obtained  by  DeLugish.  It 
is  motivated  by  the  fact  that  the  higher  radix  implementation  offers  some 
speed/hardware  trade-offs,  worth  investigating.   In  particular,  the  radix  16 
is  considered  in  four  algorithms:  division,  multiplication,  logarithm  and 
exponential.  The  central  problem  is  to  find  the  rules,  which  are  more  difficult 
when  the  higher  radix  is  used,  for  selection  of  the  digits  of  continued 
products  (sums).  The  rules  and  the  complete  algorithms  are  developed  for 
fractional  parts  in  the  conventional  range  [ l/2,  l),  of  the  floating  point 
numbers.  The  radix  l6  merely  means  "4  bits-at-a-time"  and  represents,  in 
some  sense,  the  radix  of  implementation,  not  of  an  operand.  The  exponent 
arithmetic,  being  simple,  is  not  considered.   The  use  of  a  redundant 
representation  [2,  3]  effects  the  selection  rules,  but  not  the  number  of 
steps  to  be  performed,  the  probability  of  zero,  "no  addition,"  being  too 
small  for  the  radix  l6  approach.  In  the  binary  case  [1],  the  redundancy 
is  essential  also  in  decreasing  the  average  number  of  full  steps.  As  a 
difference,  in  the  radix  l6  case,  the  number  of  steps  is  fixed  and  corresponds 
to  the  number  of  radix  16  digits,  used  to  represent  the  fractional  part.  The 
digit -by- digit  evaluation,  employed  in  the  described  algorithm,  is  not  a 
consequence  of  inherent  properties  in  the  continued  products  (sums)  approach, 
but  reflects  a  realization  strategy,  which  attempts  to  achieve  a  reasonably 
fast  implementation  for  all  functions  under  consideration,  retaining  at  the 
same  time  simplicity.  Some  comparisons  between  radix  2  and  radix  16  approaches 
are  made  and  a  more  efficient  solution,  requiring  essentially  one  "pipelined" 
arithmetic  unit  is  described. 


2.   MULTIPLICATIVE  ROEMALIZATION 

By  normalization  we  mean  a  step  by  step  transformation  of  a  given 
number  X  e[l/2,  l)  to  one  (or  to  any  other  number  N,  in  general).   For  uniformity 
and  simplicity  of  later  described  algorithms,  the  linear  convergence  is  imposed 
on  normalization.   Namely,  an  m- digit  normalized  number  is  obtained  in  (m+l) 
iterative  steps.   If  the  reciprocal  of  a  given  number  X  can  be  represented 
in  a  continued  product  form  as 

m 

1/X_  =   TT   M. 

0   i=0  x 

then  the  normalization  of  X  to  one  can  be  achieved  through  a  multiplicative 
iteration: 


\+1  =  \   '   \,  k  -  0,  ...,  m  (2.1) 


and 


m 

X  Ln  =  X_  •  TT  M.  =  1 
m+l    0   .  _  l 
1=0 


The  factors  of  a  continued  product  representation  are  of  the  form 

M=  1  +  S   .  l6"k  (2.2) 

where  S,  is  a  one-digit  constant,  so  that  an  implementation  of  (2.1)  will 
require  only  addition  and  shift  operations.   The  value  of  S  is  chosen  such 
that  the  error  e    after  step  k  becomes 


where  I S.  I    is  the  largest  constant  in  the  chosen  set  {S,  ).   In  other  words, 
1  k'max  l  k  ' 

at  every  step  k,  by  proper  choice  of  constant  S,  ,  the  k-th  digit  of  the  partially 
normalized  number  X,  will  be  forced  to  zero  or  (radix  -  l) .   Therefore,  the 
final  normalized  quantity  X    will  differ  by  at  most  16   from  unity.   To 
define  the  normalization  procedure,  one  must  find  the  rules  for  selection  of 
the  proper  value  of  S  ,  given  X,  and  the  set  [S.  }.  The  set  (S  }  is  determined 
as  follows.  To  make  the  selection  process  simple,  it  is  essential  to  have  a 
redundant  number  representation.  Therefore,  the  set  {S,  }  will  contain  more 
than  r  elements,  both  positive  and  negative.  The  maximal  value  | S,    |,  or 
better,  the  redundancy  ratio  is  obtained  from  the  practical  requirements[  3] : 
-  one  of  the  simplest  ways  to  form  the  term 
S  (X  *l6  )  in  the  recursion, 

is  to  use  a  multilevel  adder  structure,  with  corresponding 
selection-complementation  networks  generating  the 
following  sets  of  multiples: 

tO,  +1,  +2}  *  (X  'l6"k)    in  level  1; 

10,  jh,    +8}  *  (X-'l6'k)    in  level  2. 

Therefore,  the  maximal  value  should  be  I S,    I  =  10  corresponding  to  the 
redundancy  ratio  I S,    |/(r-l)  =  2/3.   The  set  of  constants  S  is  then 

{10,  ...,  I,  0,  1,  ...,  10}  (2.5) 

where  the  overbar  denotes  negative  values.  Now  the  error  (2.3)  becomes 

IVi'l  ■  ii-vMJ*f-l6"k  (2-6) 


It  is  straightforward  to  show  that  for  Xe[l/2,l),  the  reciprocal 
l/X_  can  always  he  represented  in  a  continued  product  form  using  constants 
S  from  the  set  (2.5).   Therefore,  the  normalization  procedure  defined  by  the 
recursion  (2.1)  is  always  possible. 

The  main  problem  now  is  to  define  a  practical  selection  procedure. 
The  constants  S  are  selected  in  such  a  way  that  the  error  condition  (2.6) 
is  satisfied  for  every  step.   Since  the  set  of  constants  SL  has  been  chosen 
to  be  redundant,  the  range  of  X  ,  for  all  k,  can  be  partitioned  in  the 
overlapping  intervals,  each  corresponding  to  a  particular  constant  S  .   In 
the  overlaps,  at  least  two  constants  S,  ,  differing  by  1,  are  a  valid  choice. 
S  can  be  determined  on  the  basis  of  X,  ,  but  to  keep  selection  dependent  on 
the  same  register  positions,  i.e.,  to  retain  the  same  "weights"  for  selection 
rules,  it  is  convenient  to  define  the  scaled  remainder  as 

Rk  =  l6k"1(X1^-l),    0  <  k  <  m  (2.7) 

The  selection  process  can  now  be  carried  on  the  remainders,  obtained 
recursively  from 

R^  =  161^  +  Sk  +  l6"k+1SkRk,    0  <  k  <  m        (2.8) 

This  recursion  follows  from  (2.1),  (2.2)  and  (2.7). 

If  the  general  form  of  this  recursion  for  radix  r  is  considered,  then 
the  following  remarks  concerning  implementation  requirements  can  be  made: 

i)  The  number  of  shifting  paths,  necessary  to  generate  the  last 

term  in  the  recursion  (2.8),  is  inversely  proportional  to  the 

radix  r; 
ii )  For  the  higher  radices  of  the  form 

r  =  22p,    p  =  1,  2,  3,  •  •  • 


such  that  multiple  formation  can  "be  done  with  a  cascade  of  adders, 

2 
and  with  the  set  of  S,  's  such  that  S     =  —  (r-l),  the  number  of 

k.  Kmax    j 

extra  levels  with  respect  to  the  radix  2  is  p  -  1.  Clearly,  the  full 
carry  propogation  need  be  provided  only  at  the  last  level. 
Now,  starting  from  the  condition  for  error  (2.6),  redefined  for 
scaled  remainders  as 

-  2/3  <  E^  <  2/3  (2.9) 

the  selection  rules  can  be  derived.  The  intervals  can  be  determined  by 

straightforward  calculations: 

-  for  every  S,  ,  given  bounds  on  R_    (2.9)*  find  the  interval 

boundaries  as  the  minimal  and  maximal  value  of  R.  such  that 

k 

equation  (2.8)  holds.  In  addition,  the  continuity  of  the  range 
should  be  preserved  by  retaining  only  the  overlapping  intervals. 
The  numbers  representing  boundaries  between  intervals  should  be 
simple  in  the  binary  sense  so  that  limited  precision  can  be 
used  in  implementation  of  the  selection  rules. 
From  the  definition  of  the  scaled  remainder  (2.7)  it  can  be  observed 
that  the  normalization  will  be  more  accurate  if  the  initial  step  (k=0)  is 
performed  on  X  and  then  continued  using  R.  for  k  =  1,  ...,  m.  Since 
R  =  X  -  1  this  change  in  the  procedure  is  almost  trivial. 

The  rules  of  the  selection  for  the  initial  step  can  be  made  very 
simple,  due  to  the  fact  that  S  may  be  chosen  to  be  either  0  or  1.  The 
rules  are: 

SQ  =  1   if   1/2  <  xQ  <  5/8; 

SQ  -  0    if    5/8  <  XQ  <  1. 


The  choice  of  5/8  as  the  boundary  value  is  made  so  that  the  X  will 

be  in  a  convenient  range.  From  the  rules,  the  range  of  X  is  [5/8,  5/k) . 

m  /     ,--i\ 
It  is  easy  to  see  that  even  for  m  =  1,  the  continued  products  . 7T-,  ( 1+S .  16   j 

can  represent  values  less  than  4/5  or  greater  than  8/5,  hence  making 

normalization  successful.   Therefore,  the  restriction  of  S  to  the  values 

0  or  1  is  valid. 

The  precision,  necessary  to  express  boundaries  between  intervals  is 
at  most  six  bits  after  the  radix  point.   For  convenience  the  bounds  of  the 
intervals  will  be  given  also  to  this  precision. 

For  the  step  k  =  1,  the  lower  and  upper  bounds  of  intervals  containing 

R  and  denoted  as  a  and  b,  are  calculated  for  all  possible  S  and  displayed 

in  the  following  table. 

Table  2.1 

S         a  <  64R  <  b 


10 

-26 

-23 

9 

-2k 

-22 

8 

-23 

-20 

7 

-21 

-18 

6 

-19 

-16 

5 

-17 

-Ik 

k 

-14 

-11 

3 

-12 

-8 

2 

-9 

-5 

1 

-6 

-2 

0 

-2 

3 

-1 

2 

7 

-2 

7 

12 

-3 

12 

18 

-k 

18 

2k 

-5  26  32 

-6  35  k2 

-7  k6  54 

-8  60  69 

-9  77  88 

-10  100  113 


The  values  below  the  starred  line  may  not  be  used  since  corresponding 
intervals  are  not  contiguous.  Since  -3/8  <  R-,  <  l/^,  as  follows  from  the 
initialization  rules  (k=0),  the  constants  S  can  be,  without  problems, 
restricted  to  the  set  {3,  . ..,  9)» 

By  the  same  procedure  the  interval  bounds  a  and  b  for  R  are 
calculated.  The  correspondence  between  values  of  S  and  allowed  intervals 
is  given  in  Table  2.2. 


Table  2.2 

So 

a  <  6iHL 

<  b 

2 

-         2 

10 

■42 

-36 

9 

-37 

-33 

8 

-33 

-29 

7 

-29 

-25 

6 

-25 

-21 

5 

-22 

-18 

k 

-18 

-Ik 

3 

-Ik 

-10 

2 

-10 

-6 

1 

-6 

-2 

0 

-2 

3 

-1 

2 

6 

-2 

6 

10 

-3 

10 

1*. 

-k 

11+ 

18 

-5 

18 

23 

-6 

23 

27 

-7 

27 

31 

-8 

31 

35 

-9 

35 

39 

10 

39 

k2 

The  intervals  in  which  R  may  be  found  are  contiguous  and  S  can 
have  all  values  from  the  set  {10,  ...,  10}. 

For  the  remaining  steps,  k  >  3>  the  simple  relationship  holds  between 
the  value  of  S  and  the  bounds  of  the  corresponding  interval: 


(-2Sk-l)  <  321^  <  (-2Sk+l)  (2.10) 

and 

S  e  [10,  . ..,  10} 

The  result  (2.10)  indicates,  first,  that  the  selection  rules  for  k  >  2  are 
invariant  and,  second,  that  the  selection  can  be  performed  by  rounding  the 
scaled  remainder  to  one  non-sign  digit  (in  radix  16).   The  selection  process 

itself  becomes  very  simple  after  the  first  three  steps,  due  to  the  following 

-k 
fact.   The  last  term  in  the  remainder  recursion,  l6  SVR  ,  cannot  affect  the 

most  significant  bits  of  R.  ,,,  used  in  selecting  S..  ,..,  for  k  >  2.   At  least 
°  k+1'  k+1' 

k  -  3  most  significant  digits  of  R,  remain  unchanged  except  for  possible 
complementation,  due  to  the  change  of  the  sign.   Therefore,  at  the  k-th 
step,  constants  S  ,  S   ,  ...,  S     are  known.   For  m  digit  precision,  when 

K    K+J-  ciJK— 3 

k  >  (m+3)/2  all  remaining  constants  S  ,  ...,  S  are  actually  available  and 
the  process  of  normalization  can  be  simplified.   Namely,  the  basic  remainder 
recursion  (2.8)  can  be  replaced  by  a  simple  form: 

R^  =  16R^  +  Sk,     k  >  (m+3)/2  (2.11) 

The  aforementioned  features  reveal  the  amenability  of  the  normalization 
procedure  to  higher  radix  implementations.   Once  the  process  gets  started, 
remaining  steps  are  progressively  easier. 

As  mentioned  before,  for  k  >  3  the  selection  of  S  ' s  can  be  performed 
through  rounding,  i.e.,  the  most  significant  non-sign  (radix  16)  digit  of  R 
represents  the  correct  value  of  S,  ,  after  rounding.   It  is,  then,  natural  to 
specify  the  selection  rules  for  "irregular"  steps  k  =  0,  1,  and  2  through  a 
modified  rounding  procedure  rather  than  using  a  table  look-up  or  a  direct 
combinational  approach. 


10 


The  following  definitions  are  relevant  to  the  description  of  the 
selection  rules  as  well  as  the  algorithms. 


Sign  and  magnitude  representation  of  the  constants  S,  : 

3    ± 

S.  =  (l-2s.  )  Z  s.2  ,    s.  e  [0,1}   for  all  i; 
4  i=Q  l        l 


(2.12) 


Two's  complement  representation  of  scaled  remainders: 


tan 


-i 


R,    =   -r~  +  E     r.2     ,  r.   e    [0,1}       for  all  i: 

i=l 

(m  is  the  number  of  radix  l6  digits) 


(2.13) 


Truncated  scaled  remainder: 

6    „-i 

r.2 

i=l 


\  =  "ro  +A  ri2' 


Won 


-sign  part  of  R.  : 


(2.1*0 


Then 


Tk  = 


I  r2"i 
i=l  i 


0 


6  -  0-i     •  x. 
Z     r.2       if  rQ 

i=l 


Li  t\,  =  0; 


Step-dependent  rounding  constant: 


-i 


U,  =  S  u.2  ,    u.  e  10,1}    and 

K   .  -.   1  1 

1=1 


u.  =  r.(\) 


S.  =  I  (T.  +U.  ).l6j    and 
k   L  k  k    J 


(2.15) 


(2.16) 


Sign  (Sfc)  =  Sign  (R^.) 


(2.17) 


where  [Y\   denotes  largest  integer  not  larger  than  Y. 


n 


Step  k  =  0 

For  the  initial  step,  a  modified  procedure  in  selection  is  applied. 
The  extension  of  normalization  to  the  negative  values  of  X  is  easily  achieved 
through  the  initial  step.   The  selection  rules  and  the  starting  value  of  the 
scaled  remainder  are  specified  as  follows: 
For  XQ  e  [1/2,1), 

SQ  =  1    if    1/2  <  XQ  <  5/8; 

SQ  =  0    if    5/8  <  XQ  <  1;  (2.18) 

For  negative  values  of  X  ,  it  is  a  simple  approach  to  determine  S 
by  the  rules  analogous  to  (2.18)  and  then  generate  the  negative  of  X..  "while 
calculating  R  as 


Rx  =  -3^  -  1  =  (XQ  +  XQS0)  +  2  "1U  -  1  (2.19) 

and  proceed  with  normalization  as  in  the  case  of  positive  X  „ 

Step  k  =  1 

From  the  Table  2.1,  the  interval  break  points  are  chosen  so  that  the 
simple  values  of  U  are  sufficient  to  obtain  S  .  According  to  the  proposed 
approach 


si  =  L(T1+u1)-l6J   and   Sien  (sx)  =  Sien  (\) 

In  Table  2.3  the  correspondence  between  intervals  and  rounding 
constants  (U  )  for  each  allowed  value  of  S  is  given. 


Table  2.3 


12 


\  <  ° 

sl 

6Mr1+i) 

64TX 

64U 

9 

4o 
in 

23 

22 

14 

8 

42 
45 

21 
20 

7 

44 

45 

19 
18 

10 

6 

46 
47 

17 
16 

5 

48 
49 

15 
14 

6 

4 

50 
51 

13 

12 

3 

52 
53 
54 

11 

10 

9 

3 

2 

55 
56 

8 
7 

2 

57 

6 

1 

58 
59 
60 
61 

5 
4 

3 
2 

0 

62 

63 

0 

i 

1 

0 
0 

1 

-1 

2 

3 
4 

2 
3 
4 

1 

5 

5 

6 

6 

-2 

7 

8 

7 

8 

0 

9 

9 

10 

10 

11 

11 

-3 

12 
13 
14 

15 

12 
13 
14 

15 
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Now,    it  is  a  simple  task  to  find  relations  "between  IL    and  R.    in  the 
'  Ik 

form  of  Boolean  equations.      The  derivation  of  those  equations,    listed  below, 
is  given  in  Appendix  A. 

II    =  Z     u.2"1 
i=l 

Ux  =   u2  =   0 
5  °  2  (2.20) 

U5  =  r0  +  *fk 

U6  =  *?k 

Stev  k  =  2 

Using  the  data  from  Table  2.2,  the  interval  break  points  are  selected 
in  such  a  way  that  the  corresponding  values  of  T  and  U  will  produce  correct 
values  of  S   : 


S2  =   L(T2+U2)l6j  and         Sign   (S2)   =   Sign   (Rg) 

The  correspondence  between  R~  and  U     is   shown  in  Table  2.4. 

The  Boolean  equations  are   derived  in  Appendix    B    and  listed  below: 


U.    =  U     =   u_   =   u.     =    0 
12  3  4 


u5  =  r0  +  ri(r2+V  +  r6  (2*21) 


u6  =  rQ   (r1+r2r3) 


Table  2.1* 


Ik 


\<° 

S2 

6MR2+D 

6kT2 

6i*u2 

10 

23 

1*0 

3 

26 

37 

9 

27 

36 

50 

55 

8 

31 

32 

31* 

29 

7 

35 

28 

58 

25 

6 

39 

21* 

14-2 

21 

5 

1*3 
1*5 

20 

19 

18 

2 

1* 

1*6 

17 

1*9 

11* 

3 

50 

13 

53 

10 

2 

5** 

9 

57 

6 

1 

58 

5 

- 

6l 

2 

0 

62 
63 

1 
0 

Table  2.k  Continued 
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R2  >  ° 


6kt 


6kT, 


6ifU 


0 

0 

l 

0 
1 

2 

-1 

2 

2 

5 

5 

-2 

6 

6 

9 

9 

-3 

10 

10 

r? 

13 

Jf 

ik 

Ik 

17 

17 

-5 

18 

18 

21 

21 

-6 

22 
23 
2lt 

22 

23 
24 

1* 

26 

26 

-7 

27 

27 

30 

30 

-8 

31 

31 

....    34 

^ 

-9 

35 

55 

- 

58 

38 

-10 

39 

39 

lt2 

42 

*   -  or  6kU     =  2  whenever  t/-  =  r^  =   1. 
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Step  k  >  3 

As  mentioned  before,    the   selection  procedure   for  all  remaining  steps 
consists   of  rounding: 


Sk  =  L(Tk+V,l6j>  Sign   (sk)   =  sign   (\) 


where 


Uk  =  Z  V"1  =   1^32f      i'e'> 
i=l 


u.  =  0, 


t  5 


(2.22) 


«      -1 


We  now  summarize  the  multiplicative  normalization  of  a  given  number 

!lj21  -j_ 

X     =   -x     +  Z     x.2       in  the   form  of  the  following  algorithm. 
i=l     1 

Algorithm  N   (Multiplicative  Normalization):  (2.23) 


Step  Nl.      [Initialize]  k  <-  0; 


SQ  «-  1    if  1/2  <  XQ  <  5/8; 


SQ  -  0     if  5/8  <  X    <  1; 


R1  «-XQ   (1+S0)    -  1; 


Step  N2.      [Loop] 


for  k  <  m  perform: 
k  +-  k  +  1; 


Sk  *"  l  (Tk+Uk)l6j;     Sign  Sk  «"  Slgn  Rk= 


if  k  <  (m+3)/2  then: 

-5   L    \K> 


-k+1 
\+l   -  l6\    +    Sk    +   l6  Sk\  = 


c+1 

else: 


\+i  - 16\ +  st' 


IT 


where  [Yj  denotes  the  largest  integer  not  larger  than  Y;  T  and  U  are 
defined  according  to  (2.15)  and  (2.l6)  with 

Ul  =  U2  =  °; 
u3  =  KiroV 

u^  =  \(r0+rfk)   +  Kg[r0  +  \{\   +  ??)  +  rg]  +  K 

u6  =  Ki?3r^  +  K2ro(rl+r2r5} 

and  K  ,  K  and  K  stand  for  k  =  1,  k  =2  and  k  >  3,  respectively. 

This  algorithm  will  normalize  a  given  number  X  to  one  in  (m+l) 
steps  with  the  error  bound 

|X  .  -  ll  <  2/3  •  i6"m  (2.2k) 

1  m+]_    1  _  /  -- 

and  simultaneously  generate  digits  of  the  continued  product  representation 
of  1/Xn«   The  method  of  multiplicative  normalization  is  convergent  by 
definition  of  the  procedure  and  the  existence  of  such  procedure  has  been 
shown  by  construction  of  the  selection  rules.   In  deriving  those  rules,  the 
aim  was  to  achieve  sufficient  simplicity,  not  necessarily  optimality.   The 
implementation  aspects  will  be  discussed  later.   From  the  algorithms  considered 
here,  division  and  (natural)  logarithm  are  based  upon  the  multiplicative 
normalization  and  will  be  described  in  that  order. 
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3-   DIVISION 

As  stated  before,  the  proposed  algorithms  apply  to  floating  point 
numbers  with  binary  radix  of  the  exponent.  The  radix  16  or,  in  other  words 
"k  bits  at  a  time"  is  used  to  speed  up  operations  on  fractional  parts. 
Since,  in  general,  the  exponent  manipulation  presents  no  problems,  we 
will  not  be  concerned  with  the  exponent  arithmetic  here. 

Let  Y»,  X_  e  [l/2, l)  be  fractional  parts  of  the  given  floating 
point  dividend  and  divisor,  respectively.  Then  by  multiplying  both  the 
dividend  and  the  divisor  with  the  same  sequence  of  factors  M, 

•-■J*-.-2^  (3-D 

0   Xn  f  M. 

°i=0  X 

if  XqTTM.  ->  1,  the  fractional  part  Q  of  the  quotient  is  obtained  as  Y  irM. . 
Defining  the  factors  M,  to  be  of  the  form 

\  =  1  +  Sk-l6"k  (3.2) 

and  S  e  [10,  ...,  10} 

one  can  determine  constants  S  through  the  multiplicative  normalization 
(Algorithm  N).  To  form  a  desired  quotient  Q,  let  Q.  =  Y  and  define 
recursively  the  partial  result  as 

\+l  =  \'\   =  \(1+sk*l6"k).   0  <  k  <  m         (3-3) 
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Then  Q  =  Q,  ,  has  m  correct  digits,  the  error  "bound  in  the  normalized 

divisor  being  |X  _  -  ll  <  2/3 »l6   . 

1  m+1    '  —  ' 

In  presenting  the  algorithm  for  division,  as  well  as  for  the  other 
operations,  it  is  assumed  that  the  normalization  and  the  result  evaluation 
are  carried  out  in  two  similar  arithmetic  units.   Later,  when  discussing 
implementation  aspects,  it  will  be  shown  how  the  proposed  algorithms  can 
be  realized  with  essentially  one  arithmetic  unit  with  a  tolerable  decrease 
in  performance. 

Algorithm  D  (Division) 

(AU1:  Normalization)      (AU2:  Result  Evaluation) 

Step  Dl.  [Initialize]   k  «-  0; 

Step  Nl  (Algorithm  N);      QQ  «-  Y  ; 
Step  D2.  [Loop]        for  k  <  m  perform: 

Step  IG;  \+l-\  +   Vk*"*' 

An  example  of  the  division  is  given  in  Figure  3-l»   The 

"predictability"  feature,  described  before  (2.1l),  is  apparent  at  step 

k  =  7-   the  first  five  digits  of  Rq,  when  recoded,  are  the  next  five 

constants  SQ,  ...,  S_  . 
o       12 

In  Figure  3-2,  the  basic  hardware  configuration,  consisting  of  two 
arithmetic  units,  is  shown.   The  control  part  is  not  described.  The  only 
difference  between  the  two  arithmetic  units  is  that  AU1  has  the  additional 
network  TU  and  the  five-bit  register  S,  used  in  the  selection  process. 
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k.     NATURAL  LOGARITHM 


E 
Let  X  =  X  »2   be  a  given  floating  point  number  with  fractional 

part  X  e  [1/2,1)  and  exponent  E  .   In  formulation  of  the  algorithm  in 

radix  l6,  we  follow  the  same  approach  as  for  the  radix  2  case,  given  in 

[1].   To  obtain  to  X,  the  problem  is  split  in  two  parts.   Namely, 


toX  =  toXn  +  E     to  2. 
0    x 


(h.i) 


The  algorithm  for  calculation  of  the  first  term  to  X  is  derived  from  the 


identity: 


Xn  =  X.  7T  M./  w  M. 
0    °  1*0  X  i=0  1 

-k 


ft.2) 


where  multipliers  are  M,  =  1  +  S,  »l6   and  constants  S  e  [10,  ...,  0,  ...,  10}, 
as  defined  before  for  the  multiplicative  normalization  algorithm.  Now, 


m 


m 


to  X  =  to   (X  7T  M.  )  -  to  (  $  M.) 
i=0  X       i=0  1 


(M) 


From  the  normalization  algorithm  N  (2.26)  we  know  that  the  error 
in  normalization  is  bounded  by 

|xm+1  -  l|  <  2/5-16- 


m 


Therefore,  to  (X     ir  M.  )  =  0  for  m  digits  precision  and  to  X       is  reduced  to 


^o1 


m 


m 


-i. 


toX     =   -  to  (  $  M.  )   =     S  [-  Ml+S.16      ] 
i=0  x         i=0  x 


0+A) 
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To  obtain  Ibn  X  ,  one  needs  to  perform  the  summation  of  m  +  1  sets  of 
precomputed  constants  of  the  form 

-  ton   (1  +  S.  l6"k) 

stored  in  a  fast  read-only  memory  (ROM). 

The  stored  constants  are  retrieved  from  the  ROM  using  S,  '  s,  obtained 
in  normalizing  X  .   It  is  clear  that  this  procedure  is  equally  well  applicable 
to  any  logarithm  function — just  the  set  of  stored  constants  need  to  be 
precomputed  in  the  corresponding  base.   The  summation  (k.k)   is  performed 
recursively: 

I^+1  =  1^  +  [-  en   (1  +  Sk-l6"k)]    for  k  =  0,  ...,  m 

where      L  =  0  (k.5) 

Then  L  n  =  P/n  X-.   It  should  be  noted  that  this  result  may  not  have  an 
m+1      0 

accuracy  of  m  digits  because  the  stored  constant  cannot  be  exactly 
represented  with  m  digits  and  during  m+1  additions,  errors  will  accumulate. 
If  the  result  P/r,  X  is  to  be  correct  to  m  digits  (km  bits),  i.e.,  with  error 
less  than  1/2*2   ,  then  the  precision  of  the  second  arithmetic  unit  should 
be  extended  by 

Am  =   [Pn   (m+l)/0*  21 

where  [X]  is  the  smallest  integer  not  smaller  than  X.   For  m  =  12,  this 
amounts  to  an  extension  of  k   bits.   If  the  algorithm  is  performed  in  radix  2, 
to  retain  the  same  error  bound,  the  extension  will  be  6  bits. 

The  calculation  of  the  second  term  E  Pm  2  can  be  performed  using 
conventional  multiplication  since  E  is  always  of  short  precision  compared 

-A. 
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to  X  .  The  constant  ton  2   can  be  stored  in  a  ROM.  As  an  alternative,  it 
might  be  convenient  to  implement  E  ton  2   using  a  multiplication  algorithm 
based  on  continued  sums,  as  described  in  Chapter  6.  Assuming  that  the 
length  of  the  exponent  E  is  8  bits,  such  a  solution  will  increase  the  total 
time  by  3  basic  cycles  but  reduce  the  hardware  requirements  by  an  extra  low 
precision  multiplier.  Namely,  after  2m  X  has  been  computed  in  the  second 
arithmetic  unit  (of  Figure  k.2),   this  result  is  taken  as  the  first  partial 
product.  Then  the  first  arithmetic  unit  is  initialized  to  E  ,  the  step 
counter  set  to  zero  and  the  multiplication  E  ton  2   is  performed  using  recursion 
(6.2): 

L. +1  =  1^  +  {ton  2)Sk'l6"k    for  k  =  0,  1,  2 
where 

0    m+1      0 

and  Sn,  S  ,  S  are  obtained  through  the  additive  normalization  algorithm. 
Therefore,  the  last  partial  product  P,  will  represent  the  final  result 
feX.  This  approach  has  assimilated  the  extra  add  step,  indicated  in  (^.l). 
For  higher  radices  the  set  of  constants  S  is  enlarged  and 

consequently  the  capacity  of  the  ROM  must  be  increased.  As  in  the  radix  2 

-k 
case  [1],  there  is  no  need  to  store  all  of  the  constants  [-^(l  +  S*l6  )]. 

The  possible  reduction  can  be  shown  using  the  power  series  expansion  for 

the  logarithm. 

The  constants  to  be  stored  are  of  the  form 

ton   (1  +  S  *l6"k)  =  ton   (l+a)    and 
S  €  [10,  •••,  0,  •••,  10} 
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Then 


fa  (l+a)   =  a  -  i  a2  +  H.O.T.  f  or   -   1  <  a  <  1 


2  iW^  10   -   1  +  lua  ,,,,, 

For  k>k]L,-5~8 -     ^Q+km  (k.6) 

and  m  digits  accuracy 

fa  (1  +  S.  -l6"k)   =   SL  -l6"k 

eliminating  the  need  for  storing  more  constants.   Therefore,  for  radix  16 

the  necessary  capacity  of  the  ROM  is  at  most  10m  +  15  words  of  m  digits 

(km  bits)  each,  including  the  constant  fa  2,    used  in  the  initial  step  as 

well  as  in  evaluation  of  the  term  E  fa  2. 

x 

The  algorithm  for  natural  logarithm  is  given  below: 

Algorithm  L  (Logarithm)  (^.7) 

(AU1:  Normalization)      (AU2:  Result  Evaluation) 

Step  LI.   [Initialize]    k  *-  0; 

Step  Nl.   (Algorithm  N);         LQ  -  0; 

Step  L2.   [ Loopl         for  k  <  m  perform: 

Step  N2.  if  k  <  kL  then: 

Vi~\  - "  (1  +  skl6"k)= 

else : 

Step  L3.      [Form  in  X.   +  E     fa  2] 
L  0  x  J 

Step  Al.      (Algorithm  A);  (L     =   L        ); 


0         m+1' 


-  initialize:     R  *-  E   ; 

x 


for  k  <  2  perform: 
Step  A2.  L^  -\+  (*  2)Skl6"k; 
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where  k  is  given  by  (^.6),  and  algorithm  A  for  additive  normalization  is 
given  in  (5*6).   In  the  example,  Figure  k-1,    only  the  calculation  of  the 
Sm  X  part  is  shown.   The  additional  requirement  for  implementation  of  the 
logarithm  algorithm  is  a  read-only  memory  (ROM),  connected  to  the  result 
evaluation  unit  (Figure  lj—2). 
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5.   ADDITIVE  NORMALIZATION 

Algorithms  for  multiplication  and  the  exponential,  analogous  to 
those  defined  for  division  and  the  logarithm  can  be  derived  on  the  basis  of 
continued  sums.   The  step  by  step  process,  in  which  a  given  number  X  e[l/2,  l) 


0 


is  normalized  to  zero  by  proper  choice  of  constants  S  such  that 


xn  -  §  S..16"1  =  o  (5.1) 

i=0 


where  S  e    [10,  ...,  0,  ...,  10} 


is  clearly  a  right  directed  recoding  of  X  .   It  is  termed  normalization 

as  an  additive  counterpart  of  the  previously  described  multiplicative 

normalization.   The  digits  x.  from  the  non-redundant  digit  set  [0,  ...,  15} 

are  replaced,  starting  from  the  most  significant,  with  the  digits  S 

belonging  to  the  redundant  digit  set  (10,  ...,  0,  ...,  10}.   This  recoding 

is  simple  and  exact.   Namely,  for  every  pair  of  digits  x_  x   ,  if 

x    >  10,  (maximal  allowed  value  for  S,  _)  then  this  pair  is  recoded  as 
K+l  —   '  k+l 

Sn  S.  ,n,  where  Sn  =  xn  +1  and  S,  ,.  =  -(15  -  x,  ,.).   Otherwise,  S,  =  X,  . 
k  k+l'        k    k         k+l  k+l  '   k    k 

To  define  the  selection  rules,  we  proceed  as  before.   First,  the 
scaled  remainder  is  defined  as 


Rk  =  I6k~\  (5-2) 


k-1     ,m± 
where  X,  =  X_  -  z       S.l6 

^    °   i=0   x 
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Then, 


Vi  -  l6\+i 


-k- 


^(^  -  Sk-16  *) 


l6R,    -  S  ,  0  <  k  <  m 


(5.3) 


represents  the  basic  recursion.  This  recursion  corresponds  exactly  to  the 
recursion  (2.11),  except  for  the  sign  of  S,  .  The  scaled  remainders  are 
bounded,  similarly  as  before: 


l\l  <  2/3 


(MO 


and  the  selection  rule  is  the  same  for  all  steps.  The  rule  is  simple:  S. 
equals  the  scaled  remainder  R,  ,  rounded  to  one  non-sign  digit,  the  sign 
being  that  of  R.  More  precisely,  the  selection  rule  is: 


Sk  =  ^Tk  +  V*l6j    and    Sign  Sk  =  Sign  \ 

0  <  k  <  m 


(5.5) 


where 


i=l 


-l 


if    R,  >  0; 


T  =  Z   r  2"1 
k   i=l  X 


if 


\ 


<   0     and 


\  m  "ro  +  Xri2"ls 

1=1 


U  =  1/32. 


Thi 


s  choice  of  conventional  rounding,  i.e.,  U  =  1/32  will  actually 


restrict  the  set  of  constants  S,  to  [8,  ...,  0,  ...,  8}.   This  choice  is 
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preferred  because  it  is  simplest  and  the  restricted  set  of  S  is  sufficient 
for  recoding.   The  choice  of  U  =  6/256  would  require  a  full  set  of  constants 
S.  .  Although  the  selection  rule  (5*5)  holds  for  the  initial  step  as  well, 
it  is  more  convenient  to  use  the  fact  that  for  |X  |  >  l/2,  as  is  the  case, 
|S  I  can  always  be  taken  to  be  1.  From  the  definition  of  the  scaled 
remainder  (5*2)  it  follows  that  it  is  preferable  to  start  always  with 
|S  I  =  1.   Otherwise,  for  k  =  0,  to  prevent  loss  of  accuracy,  an  extension 
of  one  radix-l6  digit  would  be  necessary  as  well  as  an  additional  right 
shift  path,  needed  only  in  this  step.  As  indicated  before,  the  result  of 
additive  normalization  is  always  exact,  i.e., 

x  -  §  s.i6_i  =  0 
i=0  x 

For  reference  purposes,  the  additive  normalization  is  summarized  in  the 

form  of  algorithm  A,  given  below. 

Algorithm  A  (Additive  Normalization):  (5.6) 

Step  Al.   [Initialize]      k  *-  0; 

s0-l; 

Ri  *~  xo  "  so; 

Step  A2.   [Loop]  for  k  <  m  perform: 

k  -  k  +  1: 


Sk  «-  l  (Tk  +  Uk)l6l;  Sign  Sk  -  Sign  J^ 

Vi  -  l6Rk  "  V 


where  U  =  1/32, 
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6.   MULTIPLICATION 

The  algorithm  for  multiplication,  given  here,  is  based  on  a 

conventional  procedure  applied  to  the  radix  l6  case.  Let  the  multiplicand 

and  the  multiplier  "be  floating  point  numbers,  satisfying  the  usual 

requirements,  i.e., 

~F 
Y=  Yn2  Y      Yn  6  [1/2,1) 


"0  0 

"cr      "0 


:  Xn2  X      Xn  €  [1/2,1) 


Again,  only  the  multiplication  of  fractional  parts  is  described, 
omitting  straightforward  exponent  arithmetic  as  well  as  postnormalization 
of  the  result. 

Consider 

P  =  Yoxo 

=  Y^  [X^  -  2  Z.  +  v  Z.l 

00  i?o  x  i=o  lJ 

where  m  is  the  number  of  radix  16  digits  in  the  fractional  parts.   If  terms 

-k 
Z  =  S  *l6   are  properly  chosen,  then 

m 

xo  -  ifo  si-16"1  ■  ° 

and 

m 

p  =  Yo  A  si'16"1  {6-1] 
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where  the  constants  S,,  as  before,  are  from  the  set  [10,  . ..,  0,  . ..,  10}. 

Applying  additive  normalization  on  X  ,  the  constants  S  can  be  obtained, 
as  described  by  Algorithm  A  (5.6). 

The  stimulation  of  the  partial  products  is  performed  simultaneously 
in  the  second  arithmetic  unit,  using  the  following  recursion: 

P.  ^  =  P.  +  Y.'S1  *l6'k,     0  <  k  <  m  (6.2) 

k+1    k    0  k     '  —   — 

where  P  -  0. 

Since  the  normalization  of  X  is  exact,  the  only  error  in  multipli- 
cation comes  from  the  single  precision  result  representation. 

The  algorithm  for  multiplication,  compatible  with  other  proposed 
algorithms,  is  given  below. 

Algorithm  M  (Multiplication):  (6.3) 

(AU1:  Normalization)     (AU2:   Result  Evaluation) 

Step  Ml.   [Initialize]     k  ^  0; 

Step  Al   (Algorithm  A);        P  «-  0; 

Step  M2.   [Loop]  for  k  <  m  perform: 


Step  A2;  P.^,  -  Pv  +  YnS.  if 


k+l  '  "k  '  "(Tk'"" 


An  example  is  shown  in  Figure  6-1.   For  implementation,  shown  in 
Figure  6-2,  an  extra  register  to  hold  the  multiplicand  is  needed. 
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7  •   EXPONENTIAL 


The  following  manipulation,  as  defined  for  radix  2  case  in  [1], 

applies  without  change  to  the  radix  l6,  producing  a  convenient  form  of  the 

v 

exponential  e   .      In  the  identity 

ex  =  gX-^e-  ^2  (7.1) 

let 

Xfo^e  =  I  +  F 

where  I  and  F  denote  integer  and  fractional  part,    respectively. 

Now 

X  Ife2Ffe2 

e     =  e  e 

JFfe2  (7.2) 

=  2   e  ' 


and  defining  X     as 


XQ  =  F^2  (7.3) 


we  obtain  the  result  in  the  form 

Y  =  YQ2  Y  =  eX  =  2Ie  °  (7-^) 

Y 

Therefore  the  problem  of  finding  e  is  substantially  reduced  to  the  problem 

X0  I 

of  finding  e  ,  the  factor  2  being  easily  incorporated  into  the  exponent 

part  of  the  result  E  .  The  exponent  X  is  any  number  such  that  Y  is  in  the 
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x 

machine  range.      Therefore,    |f|    <  1,     |x|    <to2  yielding  e       e    (l/2,    2).      It 
is  a  simple  detail  of  an  actual  design  to  obtain  F  hounded  between  -1  and 
0,    giving  e       e    (1/2,    1],    as  usual.      In  the   following  discussion,    we  assume 
that 


-1  <  F  <  0         or 

-  to2  <XQ<0  (7.5) 

The  described  approach  requires  two  extra  multiplications  before  the  main 
algorithm  can  begin.   Namely,  one  multiplication  is  necessary  to  determine 

the  terms  I  and  F  and  another  to  obtain  X  . 

X0 
The  algorithm  to  evaluate  e  ,  described  here,  is  similar  to  the 

other  algorithms,  both  in  derivation  and  in  structure. 

We  start  with  the  identity 

X^    Xn  -  0?i   (  tt  M.  )  ■+  to,   (  $  M.  )  (7.6) 

0     0      .  _  i       .  _  i 
e   =  e        1=0         1=0 

-k 

where  M=  1  +  S  •  16       0  <  k  <  m 

and  S  e  110,  ...,  0,  ...,  10} 

Once  again,  if  constants  S  are  selected  properly,  then 

X n  -M  ft  (1  +  S.16*1)]  =  0  (7-7) 

i=0      x 

and  the  result  is  obtained  as 

x      to   (  Tr  M  ) 

e  °  =  e    i=0  "   =  9   (1  +  S.16'1)  (7.8) 

i=0      x 
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i.e.,  in  the  form  of  a  continued  product.   To  define  the  selection  procedure, 

we  note  first  that 

k-1 
X.    =  X_   -     z    fc  (1  +  S.16"1),  0<k<m  (7.9) 

*         °       i=0  1 

Then 

\+l  =  Xfc  "  ^  (1  +  Skl6"K),  0  <  k  <  m  (7.10) 

and  the  scaled  remainders,  upon  which  the  selection  is  performed,  can  be 
defined  as 


E^.  =  l^"1^,  0  <  k  <  m        (7-11) 


The  basic  recursion  is 


R^       =  I6R,    -  l6k  &z  (1  +  Skl6"k),    0  <  k  <  m  (7.12) 

This  recursion  shows  that  again  precomputed  constants  of  the  form 
Qm   (l  +  S  16  )  are  necessary.  Since  these  constants  are  the  same  ones  used 
in  the  logarithm  evaluation,  all  remarks  about  storage  requirements  and 
simplification  apply  here: 

-  for  k  >  ^  «  3-\+k\  (k.6) 

the  logarithmic  constants  can  be  replaced  with  S  *l6   reducing  the  basic 
recursion  (7.12 )  to: 

R,  _  =  l6R  -  S  ,    f or  k  >  k  (7-13) 

k+1     k    k  -  1 

The  last  expression  shows  that  for  k  >  k,  the  selection  process 
will  be  identical  to  one  defined  by  the  additive  normalization  (5*6).  It 
would  be,  therefore,  natural  to  try  to  find  selection  rules  such  that  the 
similarity  with  additive  normalization  can  also  be  satisfied  for  k  <  k.. . 
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We  now  consider  rules  for  k  <  k  in  reverse  order.   First  we  recall 
that  the  selection  rules  are  determined  by  choosing  appropriate  boundaries 
between  intervals  in  R  corresponding  to  particular  S  .   Next  we  assume  that 
the  five  bit  precision  is  sufficient,  i.e.,  the  boundaries  between  intervals 
can  be  represented  as  L/32,  L  being  an  integer.   To  find  an  interval  in  R, 
corresponding  to  S  ,  the  bounds  are  determined  as 

K. 

R^  e((T^+1-l6~X   +  l6k"Vl+Skl6~k),  r^-16"1  +  l6k"^(l+Skl6"k))   (7. 1*0 

where  R  , -,  and  R    are  minimal  and  maximal  allowed  values  of  R  _ , 
respectively.   From  the  power  series  expansion,  we  have 

fcz  (1  +  S.  l6"k)  =  s.  l6~k  +  e  (7. 15) 

K.  K.  K. 


where 


1  s„2  +  1  V       1  V  +  ,_  .,, 


the   condition  of  expansion  clearly  being  satisfied. 
Let 


ek  =   l6k-ek  (7.17) 


then,  to  select  S  ,  R  must  be  in  the  interval 


or,  using  the  assumption  of  5  bit  precision 


\  «  <i  ,   J)  (7-19) 


ko 


where 

•k  ■  r2(Sk+i  +  sk  +  V1 

and  for  all  S  (7.20 ) 

\  •  ^2(\+i +  sk +  ek>J 

where  [xl  denotes  the  smallest  integer  not  smaller  than  x,  and  [xj  denotes 
the  largest  integer  not  larger  than  x.  Since  R  and  R,  have  asymptotic 
limits  -  2/3  and  2/3,  respectively,  it  follows  that  if  |  e,  |  <  1/6,  then 
a,  and  b  can  "be  determined  as 

\  -  rs(Sk+i  +  sk)i 

and  for  all  S  (7.21) 

bk  ■  i2(iW +  sk^ 

Clearly,  |e,  |  <  1/6  is  a  sufficient  but  not  a  necessary  condition. 

If  the  last  expressions  for  a,  and  b,  are  valid  then  for  selection 
of  S,  the  rounding  of  R,  to  the  most  significant  non-sign  (radix  16)  digit 
suffices,  i.e.,  the  selection  process  becomes  the  same  as  in  additive 
normalization.  Of  course,  once  S  is  obtained,  the  next  remainder  is 
calculated  using  all  terms  in  the  basic  recursion  (7-13).  By  calculating 
e,  it  turns  out  that  for  k  >  3,  |  e,  |  <  l/6  and  hence  we  have  simple  selection 
rules  as  before. 

For  k  =  2,  it  can  be  shown  that  it  is  possible  to  find  intervals 
in  R2  for  all  Sg  except  S2  =  10,  such  that  if  Rg  e  [ (2S  -l)/32,  (2S  +l)/32) 
then  S  is  the  correct  constant.  Therefore,  if  the  range  of  R  is  restricted 
so  that  S  =  10  is  excluded,  again  rounding  can  be  used  as  a  selection  rule. 


It  has  been  found  that  this  restriction  in  the  possible  range  of  R  does 
neither  affect  the  selection  process  for  k  =  1  and  k  =  0  nor  the  representa- 
tion  of  e 

For  k  =  1  intervals  are  determined,  as  before,  using  (7*1^)  and  the 
results  are  given  in  Table  7»1« 


Table  7.1 

S  a  <  32RX  <  b1 

10  15  32R 

9  lit-  15 

8  12  Ik 

7  11  12 

6  9  11 

5  8  10 

k  6  8 

3  5  6 

2  3  5 

11  3 

0-1  1 

-1  -3  -1 

-2  -5  (-11/2)  -3 

-3  32^  -11/2 

For  all  other  values  of  S  ,  i.e.,  for  (10,  . ..,  k)    intervals  are  not 
contiguous  and  hence  those  constants  may  not  be  used.   In  fact,  if  the 
possible  range  of  Kp  is  not  restricted,  S  =  k   can  be  included  in  the  set 
of  allowed  constants  in  step  1.   The  actual  selection  rules  for  k  =  1  can 
be  specified  as  in  multiplicative  normalization,  i.e.,  using  conventions 
described  by  expressions  (2.12  -  18)  one  can  determine  an  additive  constant 
U  and  through  modified  rounding  obtain  S  .  Another  choice,  which  is  given 
here,  is  to  restrict  the  range  of  R.  so  that  conventional  rounding  applies. 
This  restriction  should  not  affect  the  possibility  of  representation  of  the 
required  result.  From  Table  7.1,  if  -,171  <  R,  <  ,212  then  S  €{2,1,0,1,2,5) 


k2 


can  be  selected  applying  rounding  to  one  non-sign  digit  precision  and  no 
special  rules,  differing  from  those  for  k  >  1,  are  necessary. 

To  obtain  R  in  the  desired  range,  the  'following  initialization 
(step  k  =  0)  can  be  devised.  Since  X  e  (-  fo  2,  0]  by  assumption  and 
R  =  X  =  X  -  ton  M  the  rules  are: 


- 

Table 

7- 

,2 

xo 

Mo 

ton  M 

Rl 

[-1/8,  0] 

1 

0 

[-1/8,   0] 

[-3/8,   -1/8) 

-iA 

e     ' 

-iA 

[-1/8,   1/8) 

(-*i2,   -3/8) 

e-17/32 

-17/32 

[-.162,    .157) 

xo 

Since  e   e  (l/2,  1],  it  can  be  easily  shown  that  such  a  choice  for 
the  initialization  as  well  as  the  restricted  sets  of  constants  in  steps  1 
and  2,  i.e.,  S,  e  [2,  «>..,  3)  and  S  e  (9,  ...,  10}  can  give  the  correct 
continued  product  representation  of  the  result. 

A  summary  of  the  procedure  for  evaluation  of  e  follows: 

Preparatory  Operations  P  (7*22) 

Step  PI.        N  «-  Xfa^   e; 
Step  P2.        if  X  >  0  then: 

I  4-  [H]  +  1; 

else : 

I  -[HI? 
Step  P3 .        F  «-  N  -  I ; 

X  <-  F  ton  2 ; 


h5 

where  [N]  denotes  the  integer  part  of  N;  after  preparatory  operations,  X 
will  "be  in  the  range  (-  in  2,  0],  with  corrected  integer  part  I. 

Algorithm  E   (Exponential)  (7.23) 

(AU1:  Normalization)     (AU2 :  Result  Evaluation) 
Step  El.   [Initialize]     k  «-  0; 

Ei  -  xo  -  in  %'  Ei  *"  Mo' 

Step  E2.      [Loop]  for  k  <  m  perform: 

k  <-  k  +  1; 

if  k  <  k,  then:  K  ,,  «-  K  +  E.  S,  17~k; 

1  k+1    k    k  k    ' 

sk-^(Vuk)l6j; 

Sign  Sk  ♦.  Sign  R^; 

EL    -  l6Rk  -  l6k  Sm.   (1  -+  Skl6"k); 

else : , 

Step  A2.   (Algorithm  A); 

where  k  is  defined  in  (4.6),  and  U  =  1/32.  An  example  is  shown  in  Figure 
7-1.  The  read-only  memory,  for  this  algorithm,  communicates  with  the 
normalization  unit,  Figure  7-2.   The  remaining  configuration  is  the  same 
as  before. 
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8.   IMPLEMETWATION 

One  of  the  basic  features,  relevant  for  an  efficient  realization  of 
the  previously  described  algorithms,  is  that  the  original  operation  is 
replaced  by  the  two  much  simpler  processes  of  limited  dependency.  Through 
one  process,  the  normalization,  a  sequence  of  constants  S,  ,  the  digits  of 
the  continued  products  (sums)  are  generated.  Another  process,  the  result 
evaluation,  produces  the  final  result  using  constants  S,  .  Both  processes 
are  defined  recursively,  requiring  only  simple  hardware  operations: 
addition,  shift  and  multiple  formation. 

A  realization,  providing  one  separate  arithmetic  unit  for  each 
process,  clearly  offers  the  fastest  solution  and  the  simplest  control 
requirements.  Since  both  units  are  identical,  as  far  as  the  main  configu- 
ration is  considered,  a  replication  using  an  advanced  technology  should 
make  this  cost  acceptable.  If  the  speed  is  not  of  primary  importance,  one 
arithmetic  unit  can  be  used  in  both  processes,  performed  in  series.  For 
simplicity,  the  processes  should  alternately  use  the  arithmetic  unit  so 
that  the  current  value  of  S,  need  only  be  available.  A  "pipelining"  of 
processes  through  one  adder  and  shifting  network,  described  at  the  end  of 
this  chapter,  can  achieve  only  15^-25^  slower  operation  than  the  double 
arithmetic  unit  realization,  using  essentially  one  arithmetic  unit. 

As  mentioned  in  the  Introduction,  this  investigation  of  the  use 
of  radix  l6  in  implementation  of  the  described  algorithms,  has  been  motivated 
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by  a  possible  speed  improvement  over  the  radix  2  approach  and  by  some  trade- 
offs in  hardware  requirements.   In  the  following  comparisons  of  the  radix  2 
[1]  and  radix  16  solutions,  a  two  arithmetic  unit  realization  is  assumed  in 
both  cases.   Furthermore,  no  actual  design  being  done  in  either  approach, 
given  comparisons  are  approximate  in  nature  and  restricted  only  to  the 
size-dominant  parts.   Control  is  assumed  to  be  synchronous,  each  recursive 
step  being  performed  in  one  basic  cycle.   The  main  parts  of  the  arithmetic 
units,  used  in  comparisons,  are  as  follows. 

a)  The  adder  structure  with  the  multiple  formation  networks, 
used  in  the  radix  l6  case,  is  estimated  to  be  twice  as  complex 
as  the  corresponding  part  in  the  radix  2  case.   Namely,  in  the 
former  case,  two  adders  and  two  (l  out  2)  select  and  complement 
networks  are  required,  while  the  later  case  requires  one  adder 
with  one  select -complement  network.   The  speed  of  addition  in  the 
radix  16  case  will  be  only  slightly  decreased  if  both  adders  are 
unified  into  one  three-input  adder.   If  the  add  time  of  a  two- 
input  adder  (the  radix  2  case)  is  t  ^,  we  estimate  that  t  _/-  <  1.2t  ^, 

*  a2  a  16       a2' 

for  sufficiently  large  m. 

b)  The  shifting  network,  required  to  shift  right/left  k-digits, 
for  0  <  k  <  m  -  1,  is  simpler  for  a  higher  radix.  We  assume  that 
the  fast  shifting  network  is  realized  using  a  "barrel  switch" 
technique  [8].  Namely,  shifting  is  performed  in  two  or  more  levels 
so  that  the  combination  of  level  shifts  corresponds  to  the  required 
total  shift.   This  technique,  besides  being  fast,  ensures  low 
loading  requirements  and  the  shifting  can  be  done  using  same  paths 
both  ways:   shift  count  is  represented  in  two's  complement,  a 
negative  number  specifying  left  shift.   Implemented  in  integrated 
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circuits  technology  easily  with  its  regular  and  simple  structure, 
the  barrel  switch  as  a  standard  block  can  be  used  also  in  some  other 
operations,  e.g.,  shifting,  normalization,  etc.  Because  of  an 
additional  level,  the  shifting  network  in  the  radix  2  case  is 
estimated  to  require  30$- 50$  more  hardware  than  radix  l6  for 
m  =  48  to  6k  bits.  For  example,  if  m  =  k8,   then  level  1  may 
provide  displacements  of  0,  l6  and  32  positions,  level  2  provides 
then  displacements  of  0,  kf  8  and  12  positions,  and  in  the  radix 
2  case,  an  additional  level  3  would  be  necessary  with  the 
displacements  0,  1,  2  and  3  positions.  Speedwise,  then,  t    > 
1.3t  h-,s,   where  t  ,  denotes  shifting  delay. 

c)  The  selection  procedure  in  the  radix  2  case  requires  implementa- 
tion of  a  simple  k-   bit  comparison.  For  the  radix  l6  approach,  the 
required  precision  for  selection  is  7  bits  and  the  five  Boolean 
equations  (2.23),  costing  less  than  k-0   literals,  are  to  be 
implemented.  As  described  before,  the  selection  is  performed 
using  rounding,  so  the  additional  inputs  to  the  7  most  significant 
positions  of  the  adder  should  be  provided  as  well  as  the  5  bit 
register  S  to  store  the  current  value  of  the  constant  S,  .  Even 
with  those  requirements,  the  selection  hardware  size  is  small 
compared  with  the  rest  of  the  unit.  In  the  radix  2  case,  this 

is  even  more  true,  so  the  selection  hardware  requirements  are 
neglected  in  both  cases. 

d)  For  m  bits  precision,  the  number  of  precomputed  logarithmic 
constants,  stored  in  the  read-only  memory  (ROM)  is  about  m,  in 
the  radix  2  case,  and  about  3m,  in  the  radix  l6  case. 
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e)  The  control  part,  which  includes  also  the  step  counter  (two 
bits  shorter  in  the  radix  l6  case)  is  not  considered  as  being 
highly  dependent  on  a  particular  realization  technique. 
From  the  above  considerations,  the  hardware  requirement  ratio  for 
the  radix  2  and  the  radix  l6  is  approximately  2:3.   The  basic  cycle  can  be 
taken  to  be  the  same,  since  the  add  time  is  dominant  over  control,  selection 
and  shifting  time-   The  ROM  capacity  requirement  ratio  is  about  1:3  in 
favor  of  the  binary  case. 

Let  the  performance  of  an  implementation  in  the  radix  r  be 
P  =  Ib^r/Tr,  where  Tr  is  the  total  delay  necessary  to  evaluate  b^x   bits  of  the 
result,  as  defined  in  [k] .      Tr  is  equivalent  to  the  basic  cycle.   In  the 
radix  2  case,  the  probability  of  S,  =  0  (p  =  2/3)  is  utilized  by  providing 
an  adder  bypass  and  reducing  the  number  of  full  basic  cycles  to  m/3  on  the 
average,  where  m  is  number  of  bits.   Then,  it  can  be  taken  that  the  radix  l6 
basic  cycle  is  T  /-  =  3Tp,  since  the  number  of  basic  cycles  in  the  radix  l6 
case  is  always  the  same,  the  probability  of  S,  =  0  being  too  low.   Then, 
the  ratio  of  performances  is  P  /-/P  =  U/3  on  the  average.   If  the  efficiency 
of  the  implementation  is  defined  as  the  ratio  between  performance  and  cost 
per  bit,  then,  with  all  previous  assumptions,  E  ^-/E  ==  1  without  consider- 
ing ROM  requirements.   If  the  ROM  capacity  requirement  is  taken  into  account 
then  the  radix  2  approach  will  offer  more  efficient  design,  but  the  radix  l6 
case  will  maintain  better  performance  with  shorter  execution  time.   The 
selection  procedure  for  radix  l6  has  been  shown  to  be  sufficiently  simple. 
Even  with  the  available  efficient  technological  solutions,  the  use 
of  two  arithmetic  units  may  be  objectionable.   Since  both  the  process  of 
normalization  and  the  process  of  result  evaluation  have  addition  as  the 
basic  operation,  a  "pipelined"  use  of  the  same  adder  would  be  possible, 


50 


provided  proper  latching  of  the  operands  and  the  results  is  made.  One  way 
to  achieve  this  is  shown  in  Figure  8-1.  The  adder  with  multiple  formation 
networks  is  split  into  two  equal  parts  AS"  and  AS'  "by  breaking  the  carry 
path  and  inserting  a  one-bit  carry  register  C.  Outputs  from  the  left  half 
SN"  of  the  shifting  network  are  to  be  saved  in  a  latch  L.  The  selection 
is  carried  out  in  the  block  S  on  the  basis  of  adder  outputs  and  returns  the 
value  of  S,  .  The  initial  operands  are  in  register  B,  for  normalization, 
and  in  register  A,  for  result  evaluation.  Each  register  contains  two 
separately  controlled  halves,  B",  B1  and  A",  A'.  The  outputs  from  AS"  and 
AS1  are  connected,  under  a  separate  control,  to  the  inputs  of  A"  and  A' 
registers,  respectively.  One  separate  path  a  from  A"  to  SN"  must  be  pro- 
vided. The  operation  of  this  scheme  is  described  for  the  division  algorithm, 
with  the  help  of  Figure  8-2,  with  the  initial  control  details  omitted. 

The  normalization  process  requires  realization  of  the  recursion 

-k+1 
R    =  l6R  +  S  +  l6    S,R,  while  the  result  is  evaluated  as  Q^.+1  =  0^  + 

16"  SQ..  Since  the  operand  l6R,  in  the  first  equation  corresponds  to  the 

operand  Q.  in  the  second  equation,  additional  path  b  from  B  to  AS  must  be 

provided  as  well  as  one  k- bit  register  not  shown  in  the  scheme,  to  save 

the  most  significant  digit  of  the  left  half  of  R. 

The  operation  begins  with  the  divisor  X  in  the  B  register  and  the 

dividend  in  the  A  register.  Corresponding  to  the  scheme,  the  superscripts 

1  and  "  denote  the  right  and  the  left  half  of  each  result.  The  basic  cycle 

now  contains  two  periods,  each  period  terminated  with  the  clock  pulse. 

The  time  of  the  period  corresponds  to  approximately  one  half  of  the  full 

length  addition  time.  The  registers  are  assumed  to  be  of  master-slave  type. 

In  the  first  period,  R|  is  obtained  and  then,  simultaneously,  Q^  from  A' 

is  transferred  to  B',R'  from  AS'  to  A'  and  the  generated  carry  bit  is  saved 
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in  C.   During  this  period,  operation  on  the  left  half  is  inhibited.   In  the 

second  period,  the  right  half  evaluates  QJJ  while  the  left  half  finishes 

calculation  of  R  by  obtaining  R".   The  latch  L  and  the  path  a  provide  for 

correct  sharing  of  the  shifting  network.   Once  R"  is  obtained,  S  can  be 

determined,  the  register  transfer,  now  on  both  halves,  is  done  and  the 

process  repeats.   Shift  count  is  changed  when  Qj*  is  obtained,  matching 

requirements  of  both  recursions.   Path  b  is  selected  whenever  R  is  being 

calculated.  After  2m  +  3  periods,  where  m  is  the  number  of  radix  l6  digits, 

the  result  Q,  =  Q  ,,  is  obtained  in  the  A  register.   It  is  reasonable  to 
m+1  ° 

estimate  that  the  period  will  take  0.55-0.6  of  the  basic  cycle  so  that  the 

total  time  of  operation  in  the  pipelined  mode  will  be  increased  by  15$  to 

2.^0.      Since  the  major  blocks,  adder  with  multiple  formation  networks  and 

shifting  network,  are  reduced  from  two  to  one,  and  since  the  new  data  paths, 

latches  and  the  extra  control  are  still  significantly  less  complex  than 

the  major  blocks,  the  solution  might  be  optimal. 

Without  going  into  detailed  description,  we  mention  that  a  pipelined 

implementation  is  also  possible  in  radix  2.  With  the  assumption  that 

(t  ,.„,  +  t     ,  )  <  <  t   ,  the  total  number  of  full  cycles  will  be  about 
snxx  \j  s  s_LGC"c       ct&cL 

M/3  on  the  average,  where  M  is  the  number  of  bits,  if  bypassing  of  the 
adder  is  performed  whenever  S,  =0.   In  a  pipelined  version,  analogous  to 
the  one  previously  described  for  radix  16,  after  S  is  determined  and  before 
the  next  period  starts,  the  full  scaled  remainder  is  in  the  register  A.   If 
the  next  S  is  zero,  the  transfers  between  registers  are  inhibited  and  only 
one-bit  left  shift  on  the  register  A  is  performed,  making  selection  of  the 
next  S  possible.  When  S  becomes  non-zero,  the  normal  operation  is 
resumed.  Therefore,  the  average  number  of  full  cycles  can  be  preserved 
in  a  pipelined  version. 
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9.   CONCLUSIONS 

A  radix  16  approach  for  implementation  of  the  algorithms  based  on 
the  continued  products  (sums),  as  proposed  by  DeLugish  [1]  in  the  binary- 
case,  has  been  studied.  Those  algorithms,  in  general,  offer  simplicity  in 
mechanization  and  uniformity  in  hardware  requirements  for  a  wide  class  of 
elementary  functions.  The  use  of  a  higher  radix  makes  the  execution  of  the 
algorithms  faster,  but  additional  complexity,  both  in  the  selection  pro- 
cedures and  hardware,  must  be  considered.  For  the  radix  16  case,  it  has 
been  found  that  the  selection  rules  remain  relatively  simple.  Namely, 
the  starting  difficulties  disappear  after  the  first  three  steps,  making 
the  rules  very  simple.  Furthermore,  after  performing  the  initial  steps, 
an  increasing  number  of  constants  S  is  simultaneously  available  at  each 
successive  step.  This  property  can  be  utilized  to  simplify  normalization 
or  to  define  a  variable  radix  method,  provided  fast  and  cheap  multi -input 
adder  arrays  are  available.  Such  a  method  would  restore  fast  convergence 
of  the  algorithms,  which  is,  in  some  sense,  lost  by  specifying  algorithms 
in  a  step  by  step  mode  using  a  fixed  radix.  Hardware  requirements  for 
implementation  of  algorithms  based  on  continued  products  (sums)  are,  even 
for  the  radix  2,  greater  than  those  of  conventional  arithmetic  units  but 
still  not  prohibitive.  A  fast  variable  shift  network,  not  commonly  found 
in  conventional  arithmetic  units,  is  an  essential  part  for  the  proposed 
algorithms,  but  can  be  used  advantageously  in  many  other  operations,  like 
floating  point  normalization,  conversion  between  floating  point  and  fixed 
point  number  representations,  shifting,  etc. 
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Only  division,  multiplication,  logarithm  and  exponential  have  been 
presently  considered.  Whether  square  root,  trigonometric  and  inverse  trigo- 
nometric functions  can  be  easily  included  in  the  radix  l6  approach,  remains 
to  be  decided  by  finding  corresponding  selection  rules.   It  is  believed 
that  this  is  possible. 

Since  the  basic  hardware  is  used  for  implementation  of  many  algorithms, 
even  a  design  with  two  arithmetic  units  would  be  acceptable.   The  outlined 
"pipeline"  solution  makes  the  entire  approach  more  attractive,  since  the 
most  complex  parts,  like  the  adder  structure  with  the  multiple  formation 
network  and  shifting  network,  are  shared  by  both  processes.   This  solution 
illustrates  also  a  possible  general  approach  in  defining  arithmetic  algorithms: 
a  difficult  operation  is  decomposed  into  a  set  of  simple  processes  with  such 
interdependencies  that  the  overlapping  of  their  execution  is  feasible.   In 
this  particular  case,  the  normalization  and  the  result  evaluation  are  two 
such  processes. 
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APPENDICES 


A.   Derivation  of  equations  u.  =  f  (R,  J 


6    -i 
We  defined  U_  =  g  u.2   .  From  the  Table  2.3  it  can  be  observed  that 

i=l  X 

u.  =  f .  (r  .)    i  =  1,  . . .,  6 

J 

3  =  0,  ...,  I), 


Using  minterm  notation  we  obtain: 


-  for  r _  =  1: 


"don't  care"  minterms  are  m  ,    m  ,    ...,   m^ 


Ul  =  U2  =  ° 

u3  =  mio  +  mll  =  F2 

\  =   mi0  +  m!2  =  (F2  +  S^  (A1) 

u5  =  ro 

u6  =  m^  =  rfh 


for  r  =  0: 


"don't  care"  minterms  are  m^ ,  m  ,  ...,  ul  . 


ul  =  U2  =  u5  =  \   =  ° 

U5  =  m0  =  *fk 


IV  =   mn  =  r,r 


(A2) 


6   =   111.,   5=   i,I| 
1     3  4 
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Combining  (Al)   and  (A2),   equations  for  step  k  =  1  are: 


Ul=  U2 


u  =  r  *r 
3    0  2 


u5  =  r0  +  *?k 

B.   Derivation  of  equations  u.  =  f(R^). 

l     2 

-  for  rn  =  1: 

"don't  care"  minterm  are  m  ,  m  ,    ...,  m. 

Ul  =  U2  =  "3  =  \   =  ° 

u5  =  rQ  (Bl) 

u6  =  m5  +  m6  +  "V  +  "8  +  m9  +  m10  =  \  +  72(V^} 


-  for  r_  =  0: 

"don't  care"  minterms  are  m....,  ^o*  •••>  mir 

ul  =  U2  =  U3  =  \  =  ° 

u  =  m0  +  n^  +  m2  +  nw  +  m^  +  m  =  r^r^-Hr,)  (B2) 

u6  =  m6  +  "7  +  "8  +  m9  +  ^0  =  rl  +  r2r3 
but  u  and  U/-  can  be  given  also  as : 

u5  =  ri(?2+?3)  +  r6 
ug  =  0 

according  to  remarks  given  after  Table  2.k. 

Therefore,  for  step  k  =  2 

Ul  =  U2  =  "3  =  \  =   ° 

u^  =  rQ  +  r1(r2+r"5)  +  rg  (B3) 

U6  =  ro(rl+r2r3) 
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